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1 Introduction 

In this project a requirement for diffraction limited imaging 

for space imaging have been considered. The use of large telescopes or arrays of small 
telescopes in space are subject to time varying aberrations introduced, for example 
as a result of gravitational fluctuations. To obtain diffraction limited images from 
aberrated images requires some form of data processing algorithm to reconstruct 
the original object. 

In the past the problem of imaging astronomical objects through a turbulent 
atmosphere has been achieved using iterative algorithms on large numbers of short 
exposure, low light level data sets, to reconstruct objects from speckle type images 
[1] [2]. In contrast the space to space imaging problem requires algorithms capable 
of reconstructing objects from one or a limited number of data sets in the high level 
light regime. The proposed research was to tackle this problem in three areas. 

1. Blind deconvolution. 

2. Image reconstruction using a known aberration function. 

3. Exponential filter speckle imaging. 

The work performed upon this project has concentrated mainly upon the first 
two of the above areas. Blind deconvolution involves the retrieval of both object 
and point spread function ( PSF ) from a noisy image. This is achieved using 
iterative algorithms that use a priori information about either the object, the PSF 
or both. Over the duration of the project existing blind deconvolution algorithms 
have been implemented and new ones developed, and their relative merits accessed. 

Image reconstruction using a known aberration function has concentrated upon 
singular value analysis for the data inversion. Comparisons with many other popular 
techniques have been carried out. Testing upon real data has also been carried out 
with the use of Hubble space telescope images. 
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2 Blind deconvolution 

The intial algorithm implemented was the Aj^ers Dainty method [3] and a full testing 
of this algorithm was performed. Along similar lines a Weiner filter blind deconvo¬ 
lution algorithm and a new algorithm which used the popular Richardson and Lucy 
deconvolution technique [4,5] as its basis were implemented. Comparisons of the 
performance of these algorithms were made. 

2.1 The Ayres Dainty algorithm 

The image formed by an optical system can be represented by an integral equation 
known as a Fredholm integral equation of the first kind. 

c(a:') = [ 9 {x\x)f{x)dx (1) 

J -oc 

where 0 ( 2 ) is the image, f{x) is the object, and g{x',x) is a spatially varying point 
spread function ( PSF ). To simplify the problem the PSF is assumed to be constant 
in the image plane, i.e. isoplanatic. The above equation then becomes, 

c(x) = [ f{x')g{x - x')dx' (2) 

This form of the imaging equation allows the use of Fourier transforms so that in 
Fourier space 

C{u) = F{u)G{u) (3) 

where C{u) is the Fourier transform of c(x) and F{u) and G{u) are the Fourier 
transforms of f{x) and g{x) respectively. 

The problem of finding both object and PSF from knowledge of the convolution, 
which may also contain noise, is a difficult one. A priori knowledge of the object and 
PSF is limited, but non-negativity is an obvious constraint for incoherent imaging 
problems. Ayers and Dainty developed an algorithm to try to solve this problem 
which was essentially a generalisation of the Feinup phase retrieval algorithm [6]. 
The schematic form of the algorithm is shown in figure 1. 

2.1.1 The initial estimate of the function fo{x) 

The initial estimate of the object /o(x) is created by the use of a pseudo-random 
number generator giving positive values between zero and unity constrained to the 
extent of the convolution. 

2.1.2 Image plane constraints 

The image plane constraints upon the functions ft{x) and gi{x) help determine what 
form each function takes. A 'priori knowledge is limited, but it is known that each 
function as it is created from the inverse filter in Fourier space must be non-negative 
everywhere. Also each function must lie within the region of the convolution. The 
non-negativity constraint is implemented in the following wa}^ 

^ if /*(^) ^ ^ ^5 ( 4 ) 

" 10, otherwise. ^ ^ 
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In the Ayers Dainty algorithm the negative values of the function fi{x) are re¬ 
distributed as positive one’s to maintain energy conservation within the image space. 
In our algorithm this was found to be unecessary as Fourier plane scaling which is 
detailed later takes care of this problem. 



Figure 1. The .blind deconvolution algorithm of Ayers and Dainty 
2.1.3 Fourier plane constraints 

In each iteration the division of two Fourier transforms is required to give the next 
estimate of the function, i.e. 


F,+i(u) = 

C(u) 

G,(») 

(5) 


C(u) 

G,{u) 

(6) 


Instead of performing this operation, an averaging is used which involves incorpo¬ 
rating part of the old estimate of the function within the new estimate. Also the 
division of Fourier transforms can cause problems in areas where either jP,(u) or 
G,(u) are small. In this case inverse averaging is performed. These conditions are 
different from those given in the Ayers Dainty paper. We believe that this paper is 
in error on this point. The constraints are formalized mathematically as follows 
If |(7(u)| < (, where ^ is a noise level in the convolution, then 

= F,{u) (7) 
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If iC(u)l/l(?(u)l < F{u) then 

f,+,w = (i-mw+;3{^ (8) 

IGCu)! 

If |C'(u)|/|(j(u)| > F{u) then 

1 (l-ffl |G(u)| 

Fi{u) '^|C(u)| 

where /3 is the averaging parameter with the constraint 0 < /3 < 1. Experience has 
shown a value of ^ = 0.7 appears to give the best results. 

2,1.4 Function scaling 

It was found that a re-scaling of the functions f{x) and g{x) within the Fourier plane 
stopped a divergence in intensity of the two functions found when implementing the 
original Ayres Dainty scheme. This problem could become so bad that one of the 
functions could ‘disappear’ below the noise level. The re-scaling was performed 
before the Fourier plane constraints were applied so that 


^/m 

^>(9) 

so that Gi{u) —^ Gi{u)Sg 

(10) 

\/m 

F,{0) 

so that Fi{u) Ft{u)Sg 

(11) 


2.1.5 Results of Blind deconvolution using the Ayers Dainty algorithm 

Investigation of the effect of noise and object morphology is important for deter¬ 
mining how the algorithm will perform under different conditions. It is likely that 
objects of varying sizes'and complexity will appear in real situations, hence perfor¬ 
mance with respect to object morphology needs to be known. Also real conditions 
will produce noisy convolutions, therefore the ability of the algorithm to cope with 
noise has to be determined. Poissonian noise was added to the convolutions by tak¬ 
ing the noiseless pixel value as the mean and generating a random number obeying 
Poissonian statistics to replace this pixel value and hence create the noisy image. A 
percentage measure to give an indication of the amount of noise upon the image was 
defined as the standard deviation at the maximum pixel value in the image. This is 
then the square root of the maximum pixel value divided by that pixel value. 

As no stopping criterion exists for the algorithm an error measure, taken as the 
error between the image and the convolution of the new object and PSF, was formed. 
The best error object over a specified number of iterations was then taken as the 
result. Images are shown for the last iteration and the best error iteration. 

To test the algorithm in these conditions a simulated satellite was taken as the 
object, and a speckle type convolution generated with this object. The object and 
convolution can be seen in figure 2. 


(9) 
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Figure 2. Satellite and convolution 

In figure 3a the convolution is noiseless, and both the last iteration and best 
error objects show a reasonably clear satellite, although they appear under the PSF 
heading. This is to be expected as the algorithm makes no distinction between 
the two, also the object and PSF can swap as the algorithm iterates. Up to 150 
iterations the errors are very large, but after this point, apart from a few spikes the 
errors remain comparatively low. 

An investigation of the performance of the algorithm in the presence of noise was 
undertaken. At 0.01% noise in figure 3b the satellite appears remarkably well. ^ 
the last iteration results fringes can be seen across the images, especiafiy the PSF, 
this is characteristic of the change over between object and PSF, and is associated 
with a large error. This is difficult to see in the error graph, but looking closely 
reveals a large spike beginning to appear at the last iteration. Further increases m 
the noise did not produce recognizable objects. 

2.1.6 Conclusions ' , 

It has been shown here that it is possible to obtain some good reconstructions 
from the Ayres Dainty algorithm, but as was shown in the first quarterly report the 
speckle pattern image of the satellite was probably a special case. Also the tolerance 
to noise was especially poor, therefore other algorithms were considered in light of 
these drawbacks. 
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Figure 3a. The blind deconvolution of a satellite from a speckle pattern with zero 

noise 
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Object, and best error object 



PSF, and best error PSF 


Number of Iterations 


Figure 3b. The blind deconvolution of a satellite from a speckle pattern with 

0.01% noise 






















2.2 Blind deconvolution using a Weiner filter approach 

To improve the noise tolerance of the algorithm it was thought that a Weiner filter 
would be more appropriate than the complicated division in Fourier space method 
of the Ayres Dainty algorithm. The division in Fourier space is carried out using 
the modified formula, 

where $(u) is the optimal or Weiner filter, and is given by, 

$ (u) =_- (13) 

|5(u)P + lAr(u)p ^ ^ 

with S{u) = F{u)G{u) and where N{u) is the noise spectrum. Hence 

G-{u)C{u) 


F{u) = 




The |7V(ii)p/lF(u)p term is unknown and is replaced by a constant fi to be deter¬ 
mined by experiment. 


2.2.1 The Blind Weiner Algorithm 

A blind deconvolution algorithm is constructed in the same manner as before but 
the Fourier plane constraints are replaced by the formula 


F(u) = 


when determining F'(u), and 


G{u) = 


iG(u)P + /z 


F^{u)C{u) 

|F(u)12 + /X 


when determining G[u). The image plane constraints and the function scaling in 
Fourier space remain the same as those used in our form of the Ayres Dainty algo¬ 
rithm. This algorithm is similar to one given by Davey ei al [7], but there algorithm 
assumes further a priori knowledge in the form of a known object support. 


2.2.2 Results 

In contrast to the Ayers Dainty algorithm the performance of this algorithm upon 
speckle type images was particularly poor. It was not possible to obtain any recog¬ 
nizable reconstructions. However, convolutions that used single Gaussians for the 
PSF faired better. The original object, Gaussian PSF and convolution can be seen 
in figure 4. 

The Ayers Dainty algorithm performed badly upon this type of image, but as 
can be seen in figure 5a definite cross and Gaussian reconstructions have been found. 
It was further found that with this type of image and about 1.5% noise, i.e. ~ 4000 
photons in the maximum pixel value of the image, that reconstructions of a similar 
quality could be obtained, figure 5b. 
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b. Reconstructions with 1.5% noise. 


Figure 5. Blind deconvolution using the Weiner filter algorithm. 


2.2.3 Conclusions 

This section has shown that it is possible to achieve blind deconvolution in the 
presence of more realistic levels of noise, but for real data it is likely that noise levels 
will be greater than 1.5% will be encountered. Also the algorithm is stiU susceptible 
to the type of image given to it, hence further improvement is still needed. 

2.3 Blind deconvolution using the Richardson and Lucy 
algorithm 

The Richardson Lucy deconvolution algorithm has become very popular in the fields 
of astronomy and medical imaging. Initial!}^ it was derived from Bayes theorem in 
the early seventies by Richardson and Lucy [4,5]. In the early eighties it was red¬ 
erived by Shepp and Vardi [8] as an algorithm to solve positron emission tomography 
(PET) imaging problems where Poissonian statistics are dominant. Their method 
used a ma>dmum likelihood (ML) solution, which was found using the expectation 
maximization (EM) algorithm of Dempster ei al [9]. The reason for the popularity 
of the algorithm is its implementation of ML and its apparent ability to produce 
reconstructed images of good quality in the presence of high noise levels. We there¬ 
fore assumed that a blind form of this algorithm would have the same characteristics 
[10]. A blind deconvolution algorithm similar to the one shown here has also been 
developed by Holmes [11] using the EM algorithm of Dempster [9]. 

To begin with a brief review of the Richardson Lucy deconvolution method will be 
given, then the blind form of the algorithm will be presented. Since Bayes theorem 
relates conditional probabilities the resulting algorithm takes into account statisticiil 
fluctuations in the signal, and therefore has the ability to reconstruct noisy images. 
Bayes theorem is :- 

where P{y\x) is the conditional probability of an event y given event x. P(x) is 
the probability of an event x, and P{x\y) is the inverse conditional probability i.e. 
the probability of event x given event y. The probability P{x) can be identified as 
the object distribution /(x), the conditional probability P[y\x) can be identified as 
the PSF centered at x i.e. g{y^x) and the probability P[y) as the degraded image 
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or convolution c(y). This inverse relation allows the derivation of the iterative 
algorithm given below. 


.w = / 


g{y,x)c{y)dy 


where i is the iteration number. If an isoplanatic condition exists then this equation 
can be written in terms of convolutions:- 


fi{x)®g{x) 


where ® is the convolution operation. The PSF ^(a;) is known, so the object /(x) is 
found by iterating the above equation until convergence. An intial guess is required 
for the object fo{x) to start the algorithm. Then in subsequent iterations, due to 
the form of the algorithm, large deviations in the guess from the true object are 
lost rapidly in initial iterations whilst detail is added more slowly in subsequent 
iterations. Advantages of this algorithm include a non-negativity constraint if the 
initial guess fo{x) > 0. Also energy is conserved as the iteration proceeds, which is 
easily seen by integrating both sides of equation 18 over x. 

In the blind form of this algorithm two of these deconvolution steps are required. 
At the blind iteration it is assumed that the object is known from the — 1 
iteration. The PSF g^{x) is then calculated for a specified number of Richardson 
Lucy iterations as in equation 20, where the i index represents the Richardson Lucy 
iteration. This equation is essentially an inverse of equation 19 since the object and 
PSF have reverse roles and it calculates the PSF from the object. Then f^{x) is 
calculated for the same number of Richardson Lucy iterations. This is done using 
the PSF evaluated from the full iteration of equation 20. In this case the iteration 
is performed in the normal manner of equation ^9^,as shown in equation 21. The 
degraded image is again given as c(x) in both equations 20 and 21. The loop is 
repeated as required. 














( 21 ) 


The above equations are shown in one dimension, the extension for two dimen¬ 
sional images is straightforward. Starting guesses are made for the object fo{x) and 
the PSF algorithm loop of the form of figure 6 is performed. No posi¬ 

tivity constraints are required because the above equations always ensure this. The 
algorithm [10] is different from the Holmes [11] algorithm as only two Richardson 
Lucy iterations are performed within one blind iteration, one for an object evalua¬ 
tion, and one for the PSF evaluation. It was found that the simulated images used 
did not perform well with this type of iteration, but when the number of Richardson 
Lucy iterations with one blind iteration was increased to about ten a much better 
performance was obtained. 





Figure 6. The blind deconvolution algorithm based upon the Richardson and Lucy 

deconvolution algorithm 

2.3.1 Results 

In the case of the cross convolved with a Gaussian which gave reasonable results 
for the Weiner filter algorithm much improved results have been obtained in terms 
of noise performance. Reconstructions have been obtained when the photon noise 
was 10%, which allows the possibility of using this algorithm on real data. Ihe 
convolution can be seen in figure 7a. and the reconstructions of the cross and the 

Gaussian PSF in figure 7b. , . r -i j * 

This algorithm was also apphed to a speckle image of the cross but failed to 

produce any recognizable reconstructions. The algorithm was applied to many other 
images where gaussian PSFs were used and it was found that as long as the 
of the PSF was not to. severe then reasonable reconstructions could generally be 
obtained, and in some cases with noise levels as high as 15%. 



i. The convolution with 10% noise. 


12 






b. Reconstructions of the object and PSF. 

Figure 7. Blind deconvolution using the Richardson Lucy algorithm. 

2.3.2 Conclusions 

As hoped this algorithm has the ability to produce reconstructions at noise levels 
that are likely to be encountered in real situations. The quality of the reconstructions 
is good, but the object still significantly blurred. The way to proceed would be the 
used of further a priori information. 

3 Semi - Blind deconvolution 

It was thought that further a priori information could be incorporated by assuming 
partial knowledge of the form of the PSF. In a real situation it may be known that 
a telescope suffers from spherical aberration, but due to time varying factors such 
as a changing gravitational field the extent of this aberration may not be known. 
This would reduce the number of unknown variables in the deconvolution from 
maybe thousands of pixel values, to one or two unknown constants. We termed this 
approach Semi-Blind deconvolution [10]. 

3.1 The Weiner Semi-blind algorithm 

This algorithm used the Blind deconvolution using a Weiner filter as its basis. The 
only part of the algorithm which was altered was the image plane constraints on 
the PSF. In the blind algorithm this was just non-negativity within the convolution 
area, and zero outside it. This was replaced by a least squares fitting procedure. 
Initially convolutions were created with Gaussians, so Gaussians of var 3 dng widths 
were compared to the function input to the routine. The Gaussian giving the least 

error in fitting was then chosen as the PSF. 

To illustrate how well this algorithm performed, figure 8 shows the reconstruc¬ 
tions of the object and PSF at every iteration. 
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Figure 8a The random guess start for the object, and the convolution 



Figure 8b The results from the first iteration 



Figure 8c The results from the second iteration 



Figure 8d The results from the third iteration 
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This algorithm has converged within three iterations and has produced a perfect 
reconstruction of the satellite. The initial convolution was in fact noiseless. When 
trying to reconstruct with noisy images the algorithm always converged upon the 
delta function solution i.e. the Gaussian of smallest possible width. Therefore the 
impressive result obtained in figure 8 would never be found when using real data. 

3.2 Semi-blind deconvolution using the Richardson and 
Lucy algorithm 

To try and alleviate the noise tolerance problem the Richardson and Lucy algorithm 
was tried yet again. The semi-blind form of the algorithm took as its basis the blind 
algorithm, but in this case a number of blind iterations were performed before fitting 
a Gaussian to the function it had found for the PSF. 

The results for this algorithm have shown remarkable noise tolerance. In figure 
9 results are shown for semi-blind deconvolution on a series of point sources. The 
image contained approximately 20.0% noise. 



Figure 9a The object and image with 20.0% noise. 



Figure 9b The reconstruction of the object and the chosen fitting Gaussian PSF. 

The algorithm was tried on the noisy image of the cross used earlier for the pure 
blind deconvolution work. Although the PSF was fitted in each iteration with a 
Gaussian of the correct size the results were not very good, in fact the pure blind 
deconvolution results were better. Therefore it was decided that a Gaussian fitting 
process should be performed after the blind deconvolution, and then a specified 
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number of Richardson-Lucy iteration performed with the guessed Gaussian. The 
step width in Gaussian fitting was obviously important, with the step width been 
one pixel for the Gaussian radius at the 1/e height the correct PSF width of 3 pixels 
was guessed. At a step width of 0.1 pixels the Gaussian width obtained was 3.2 
pixels. The results for both cases can be seen below. 



Figure 10 The reconstructions of the cross of figure 4 with Gaussian step widths of 

1.0 and 0.1 pixels. 

The results show that the slight error made in finding the width of the Gaussian 
PSF do not make the reconstructions significantly worse. This is probably due to the 
high level of noise upon the image resulting in a loss of a large amount of information. 
To show the impressiveness of these results straight forward deconvolutions using 
Fourier regularisation and the Richardson Lucy are shown in figure 11. 



Figure 11 Reconstructions using (a) Fourier regularisation (b) Richardson and 

Lucy. 


It can be seen that the semi-blind deconvolution results are comparable to the 
usual methods of deconvolution. 

To extend this work to use in realistic situations more than one fitting variable 
may be needed to describe the PSF accurately. To test this a simple PSF was 
created with the functional form 


j/(^) = 


k 


.4*exp(1.0)r- 
-TTl- + tSk 




( 22 ) 
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where r is the radius and the variables Ak,Bk,Ck were given the following values:- 


= 0.0, Ai = 0.05 

= 1 . 0 , B2 = 0.0 

C, = 1.0, C 2 = 5.0 

The variables A 2 ,Ci,C 2 were allowed to change their values, so that the PSF was 
a Gaussian plus a Gaussian times its radius squared. Incorrect values for these 
variables were entered into the program and a PSF created. Then as before a blind 
deconvolution process evaluated a new object and PSF. A PSF with the above func¬ 
tional form and free variables ^2,(7], (72 was fitted to the evaluated PSF using a 
Levenberg-Marquardt non-linear, least squares fitting routine. This returned new 
values for A 2 ,C^,C 2 , and the process was repeated for a specified number of itera¬ 
tions. The PSF and noiseless convolution can be seen in figure 12, the object was 
the cross shown earlier. 



Figure 12 (a) PSF (b) Convolution. 

To compare the results of this semi-blind deconvolution algorithm a Richardson 
and Lucy deconvolution was performed using the known PSF. The result of this can 
be seen in figure 13. 



Figure 13 The Richardson Lucy deconvolution of the above image with the above 

PSF 

This can be compared with the results of 15 iterations of the semi blind decon¬ 
volution algorithm. 






Figure 14 Semi-blind deconvolution (a) The object (b) The PSF. 

The evaluated parameters for the PSF are = 0.0508, Ci = 1.05,(72 = 5.07 
with 1.7% error in the evaluation of the PSF. In the case of noisy images the results 
in terms of fitting the correct PSF were just as good. Photon noise of about 4.0% 
was added to image in figure 12b. This image along with a Richardson Lucy decon¬ 
volution are shown in figure 15. 1000 iterations were used in this reconstruction. 



Figure 15a The image with 4.0% noise 



Figure .151) The Richardson Lucy deconvolution of the above image at a 4.0%) noise 

level with 1000 iterations. 

In figure 16 the result of 15 iterations of the semi*blind deconvolution method 
are shown. The values for the fitting parameters were .4i = 0.0511, (7i = 1.056,6*2 = 
5.072 with an overall fitting error on the PSF of 1.7%. 
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Figure 16 Semi-blind deconvolution (a) The object (b) The PSF. 

In both of the above cases the semi-blind deconvolution algorithm showed definite 
signs of convergence. A blind deconvolution algorithm was run on the noiseless image 
and this did not converge to the correct result, but it did pass with 2.3% error of the 
correct PSF. The semi-blind algorithm was also run upon images containing 7.0% 
and 10% noise. Again these tests converged but to the wrong values for the fitting 
parameters, but it could be seen that at about 5 iterations the fitting parameters 
were correct. The convergence properties are shown in more detail in the next series 
of tests. 

A second PSF was chosen and an image created. In this case the fitting parame¬ 
ters were sUghtly different Aj = 0.1, = 1.0, Cs = 5.0. This made the image rather 

more blurred. 



Figure 17 (a) PSF (b) Convolution with 1.0% noise. 

In this case convergence was found for the noiseless case and the 1.0% noise 
image. The starting values entered to the program were A 2 = 0.5, Ci = 3.0, C 2 = 7.0 
giving a 74.0% error in the PSF. It can be seen in figures 20 and 21 that at this 
noise level the algorithm is converging upon the correct values. To obtain these error 
graphs the true PSF must be known which is slightly false since the true PSF will 
not be known in a real situation, but it does show the convergent properties of the 
algorithm. This does not occur in cases such as the Ayers-Dainty, blind Weiner and 
blind Richardson Lucy algorithms, see figures 3a and 3b for example. In the case of 
the blind Richardson Lucy algorithm the image in figure 17b was used to compare 
with the results of the semi-blind algorithm. It was found that the algorithm didn t 
converge to the correct values for the fitting parameters, and in fact the algorithm 
eventually diverged. 
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The final \'alues for the fitting parameters of the semi-blind algorithm were = 
0.102, Cl = 1.06, C 2 = 5.03 with an overall error in the PSF evaluation of 1.09%. 
The results for the 1.0% noise image can be seen in figure 19. Figure 18 shows the 
reconstruction using Richardson Lucy with the real PSF, the results can be seen to 
compare well with those of the semi-blind deconvolution in figure 19. 



Figure 18 The Richardson Lucy deconvolution of the above image at a 1.0% noise 
level, with the above PSF. 1000 iterations were used. 




Figure 20. Variation of the fitting parameters with iteration number for 1.0% noise. 
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At 2.0%, 3.0% and 4.0% noise levels similar results were obtained and conver¬ 
gence seen at the correct values of the fitting variables. The results of the 4.0% 
noise image are shown in figure 23. Figure 22 shows the Richardson and Lucy 
deconvolution after 1000 iterations. 



Figure 21. The error in the PSF with iteration number for 1.0% noise. 



Figure 22 The Richardson Lucy deconvolution (a) The image with 4.0% noise (b) 
The reconstruction with 1000 iterations. 






Figure 23 Semi-blind deconvolution with 4.0% noise (a) The object (b) The PSF. 







Again the results compare quiet well. Figure 24 shows the variation in the fitting 
parameters and figure 25 the percentage error in the PSF with iteration number. 
Again convergence is seen. When the same image with 6.0% noise was tried, the 
results were not so good. Figure 26 shows the variation of fitting parameter and 
figure 27 the percentage error in the PSF and it can be seen that convergence 
is reached after 8 iterations and then the algorithm starts to diverge. It would 
therefore appear that the algorithm has a certain noise tolerance. 



Figure 24. The variation of the fitting parameters with iteration number for 4.0% 

noise. 



Figure 25. The error in the PSF with iteration number for 4.0% noise. 
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Iteration number 


Figure 26. The variation of the fitting parameters with iteration number for 6.0% 

noise. 



Iteration number 


Figure 27. The error in the PSF with iteration number for 6.0% noise. 

3.3 Conclusions 

In many real situations it may be the case that some knowledge of the PSF can 
be obtained. Therefore functional forms for the PSFs were chosen with a number 
of unknown variables. It was found that accurate deconvolutions could be made, 
very near the quality of a deconvolution with full knowledge of the PSF. It appears 
that the extra a priori knowledge that the semi-blind algorithm uses stablises the 
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algorithm and gives it some convergent properties. It is hoped that this work can 
be extended to real images with PSFs containing unknown amounts of aberration 
where the algorithm will evaluate both the aberration coefficients and the object. 

4 Image reconstruction using a known aberra¬ 
tion function 

Image restoration has been a subject intensively investigated over many years and 
with the increased power of computers available today remains a topic of much 
interest. The most recent and well pubbeised example of this has been the restora¬ 
tion of images from the Hubble space telescope (HST), before the corrective optics 
were added. In the case of the HST and many other imaging problems the Point 
Spread Function (PSF) is known. To obtain high-resolution images then requires 
the solution of the linear integral equation that relates the object to the image. 
This integral equation is known as a first-kind Fredholm equation and is a classical 
example of an ill-posed problem, namely a problem whose solution is affected by 
numerical instability and therefore is strongly noise dependant [12]. Approximate 
and stable solutions can be obtained using suitable techniques. 

In the case of low-aperture systems an isoplanatic approximation can be made 
and the Fredholm equation becomes a convolution equation which is far simpler to 
solve. Many algorithms exist to solve this problem, there are the direct linear meth¬ 
ods such as Tikhonov regularisation, Landweber and singular value decomposition 
(SVD). Among the non-linear iterative techniques are the Landweber iteration with 
imposed positivity, conjugate gradients and the Richardson and Lucy method [4,5]. 
All of the above methods except SVD generally employ the Fast Fourier Transform 
(FFT) for efficient computation. This allows restoration of large images because of 
the nlog{n) speed scaling of the FFT. 

In the case of real data no simple quantitative measure can be made of algo¬ 
rithm performance. Therefore simulated HST images were used to gain an accurate 
comparison of the algorithms. This work was published at a Hubble space telescope 
workshop in November 1993 [13]. A resume of this work is given here. 

k number of methods were used to compare and contrast their performance. 
The methods considered were the following:- 

1. Tikhonov regularisation 

2. Landweber iteration 

3. Landweber iteration with imposed positivity 

4. conjugate gradients 

5. conjugate gradients with imposed positivity 

6. Richardson and Lucy 

Allofth ese methods contain a free parameter which is the so called regularisation 
parameter. In the first method this is essentially a parameter that reduces higher 
frequencies in the solution, and in the other methods it is the number of iterations. In 
the case of simulated images the free parameter can be chosen as the value that gives 
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the minimum distance in terms of euclidean norms between the approximate solution 
and the true solution. This then gives the best solution. In the case of real data 
this obviously cannot be done, so one has to choose some empirical criterion which 
depends upon the problem been investigated. For the simulated data used here 
the minimum in distance between the euclidean norm of the true and approximate 
solutions will be used to give the optimum regularisation parameter. 

The approximate solution to the problem g = Af + 77 is given by:- 

+ /xl) A^g (23) 

where is the regularisation parameter. In case 1 Fourier transforming leaves:- 

m / A ^ (^) ^04^ 




|Ff(u;)P + /x 


where F),(w) and G{u)) are the Fourier transforms of and g. H{uj) is the Fourier 
transform of the PSF. In case 2 the iteration scheme is given by:- 

fo = 0; f„+i = f„ - Tr„ (25) 

where 0 < r < 2||A||“^ is the stationary relaxation parameter and 

Tn = A^Afn - A*g (26) 

This iterative method is equivalent to filtering since the result of n iterations can be 
written explicitly in Fourier space as:- 

= [1 - (1 - .liXMP)"] PT) 

In method 3 this cannot be used since the imposition of positivity can only be carried 
out in real space, hence this algorithm requires the use of two FFTs per iteration. 
The relaxation parameter r for both methods can be set as follows:- 


where |-ff(ci;)lrnnT — Tnax.^>\H\. Finally the conjugate gradients iteration scheme is 

fo = 0; fn+i=fn“r,p„ (29) 


where 


Pn = ) 

with Po = To = -Ag, = llr„l|Vl|r,_ill' and the non-stationary relaxation 

parameter = (Th) Pn)/|| Ap„|p. This iteration scheme can easily be written in 
Fourier space. The last method Richardson and Lucy was described earlier. 

In the numerical experiments two test problems were considered. One of the 
objects w’as a simulated elliptical galaxy and the other was a star cluster. The sim¬ 
ulated blurred images were obtained by convolving the object / with a calculated 
PSF and adding Poissonian and readout noise. The discrepancy between the recon¬ 
structed solution and the true solution / may be measured by the RMS norm 
e of their difference divided by the norm of /. 

MfT zIil (31) 
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Algorithm 

iteration number 

RMS 

CPU time 

Richardson Lucy 

1006 

0.0606 

21 min 18 sec 

Landweber with pos. 

1038 

0.0484 

12 min 58 sec 

Landweber 

direct 

0.1020 

2 sec 

Conj. grad. | 

no 

0.1355 

2 min 5 sec 

Conj. grad, with pos. 

100 

0.1027 

1 min 52 sec 

Tikhonov 

direct 

0.1115 

1 sec 


Table 1: Algorithm comparison - Star cluster image 


For the direct methods the regularisation parameter is adjusted to minimize c and 
the iterative methods are run until a minimum in e is found. 

The best performance for the restoration of the star cluster was provided by the 
Richardson Lucy and the Landweber with positivity methods. They gave accuracies 
of 6% and 5% respectivel 3 \ The full results are shown in table 1. 

These results would appear mathematically to be exceUent but practically the 
results are far from satisfactory as many false stars appear when examining in detail. 
This is partially due to the fact that the solution itself is represented by discrete 
points, whereas the image and reconstructed solution are blurred representations 
thereof. In fact the choice of such solutions are particularly demanding for the 
restoration routines since they are delta-function like objects with sharp cut-offs. 
As a matter of fact the Landweber with positivity seems to perform better, this is 
highlighted by table 1 where it is shown that this method requires far less compu¬ 
tational effort than the Richardson Lucy method. The other methods seem to be 
fast but this is at the expense of accuracy. Some example restorations are shown in 
figure 28 for the Richardson Lucy and Landweber with positivity methods. 

In the case of the elliptical galaxy a similar set of results were obtained but 
this time the Richardson Lucy method seemed to out perform the Landweber with 
positivity method. 

The mathematical problem arising from HST image restoration provides an ex¬ 
ample of inverse problems which can be fruitfully treated by regularisation tech¬ 
niques. Among these methods Landweber with positivity seems to provide the best 
compromise between accuracy and computational effort. On the other hand it is not 
evident from the numerical experiments performed here which is the best choice, a 
larger set of test problems has to be taken before making firm conclusions. 

The HST was our main source of real data. One such real image was the super¬ 
nova snl987a, algorithms such as Tikhonov regularisation, Richardson Lucy etc were 
applied to this image. It was found that Tiknonov regularisation gave some very 
poor reconstructions, whereas the Richardson and Lucy method was quite impres¬ 
sive. In figure 29a the supernova image is shown on a log scale and the Richardson 
Lucy reconstruction is shown in figure 29b. 
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Figure 28 a. The simulated star cluster image. 






Figure 28c. The Richard 



Figure 28d. The Landweher 






r 



Figure 29a. The HST image of the snl987a supernova. The intensity is scaled 

logarithmically. 



Figure 29b. The Hichardson Lucy 

intensity is 


reconstruction of the snl987a 
scaled logarithmically. 


supernova. 


The 
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5 Image reconstruction using scanning methods 

A method of image reconstruction that does not rely upon the fast Fourier transform 
is SVD. To perform an SVD on an n x n matrix requires the order of operations, 
and since images require the inversion of an operator of size li^ x SVD scales as n^. 
It is therefore readily apparent that the use of SVD upon large images is completely 
out of the question not only because of the storage reqiiirements of an x matrix 
but the also the extremely large times required to perform the SVD. On the other 
hand once the singular system of the optical arrangement has been calculated and 
stored, efficient high-resolution-imaging can be obtained for that system for any 
object. Also the SVD of the imaging system is preferable since it produces the most 
natural mapping from object to image space, and hence allows a suitable method 
to reconstruct images contaminated with noise. This is achieved with the use of a 
regularisation parameter or by cutting out the singular values which are below the 
signal to noise ratio. As a result of these points we feel that SVD methods are one 
of the better ways for of restoring images. Therefore a method was required that 
allowed the use of SVD and was computationally efficient. 

5.1 Introduction 

A technique was devised that scanned the image and reconstructed in parts. This 
relied upon the local nature of the PSF, i.e. the PSF is much smaller than the 
image. In most cases this is true. 

We show that this technique is vastly superior to SVD in terms of the number 
of operations required. Also we show how the blocking effects that would normally 
be expected from the reconstruction of an image by parts are removed by the use 
of our scanning method. Reconstructions are then performed upon large complex 
images. Finally reconstructions are performed upon images formed with a spatially 
variant PSF. 


5.2 Singular value decomposition in diffraction limited 
imaging 

The linear imaging equation that we wish to solve to restore an image is the first 
kind linear Fredholm integral equation which can be written as, 

9{x) = {Hf){x) = J h{x,x')fix')dx’ (32) 

where the vectors x and are two dimensional vectors with components (xi,X 2 ) 
and {x\^X 2 ) respectively in the data and object spaces. To simplify this equation an 
approximation often made is the isoplanatic approximation, which leaves equation 
32 as a convolution. 

5(x) = {Hf){x) = I h{x- x')f{xyx' (33) 

The imaging equation can also be formulated in the discrete case as a vector-matrix 
problem. 

g = H.f (34) 


30 


The solution of the imaging problem is the inversion of the integral operator in 32 
and 33 or the matrix operator in 34. One method of approaching this problem is 
an eigenvalue or singular value decomposition. Here we will give a brief resume of 
some points relevant to image reconstruction via SVD. 

If an operator A is compact and self adjoint then it maps as A : A A and a 
set of eigenvectors {uk} form an orthonormal basis in the closure of the range of A 
i.e ^(A). Then for a vector in iZ(A), 

Ai = ^k{f,Uk)x'^k (35) 

k=0 

where (‘,*)x inner product on space X and Xk are the eigenvalues of A. If 

the operators 

A = : X ^ X 

A^ = HH^ :Y -^Y 


have common non-zero 

eigenvalues cr^ and sets of eigenvectors {ujt} re- 

spectively then 


H^Huk = (rluk 

(36) 



(37) 

where k = — 

1. Therefore we can always choose eigenvectors Uk 

and Vk 

such that 


Hufr = (XkVk 

(38) 


= cTkUk 

(39) 

Hence we have 


Hf = ak{f,Uk)xVk 

(40) 


kzzO 

ii^g=Y,^k{g,Vk)yUk 

k=0 

(41) 


The system is called the singular system. The singularvectors {ujt} span 

the object space and the singularvectors span the image space and the mapping 
from one to the other is controlled through the singular values (Tk* A solution still 
requires the inversion of the operator H. An approximate solution can be found as 
the unique least squares solution of minimal norm (Moore - Penrose solution) i.e, 


k=0 

The propagation of relative errors from data to solution is controlled by the condition 
number. If is a small variation of g the corresponding variation of f is then, 

- < cond[A)- j~— (43) 

V ilgllr 



where cond{A) = aojoTR^i, If cond[A) is large the problem is said to be ill- 
conditioned and small variations of the data produce a completely different solution 
i.e. the solution is unstable to noise. 
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In order to find an estimate of the unknown object the effect of noise must be 
taken into account. Hence the true imaging equation is 

g = Hf + 7? (44) 

where rj £ Y is the noise vector. A signal to noise ratio can then be defined as 



(45) 


This requires that the smallest singular value of H, be greater than a. If this 

condition is not satisfied then only those singular \'alues that do satisfy this condition 
can be used in the reconstruction, i.e. > a and the approximate solution is given 

by ^ ^ 

f H — (g.i'AOvWjt (46) 

k=0 

In fact this is a linear filter to the Moore-Penrose solution 

f = (47) 

A-=() 

where the filter contains a step function that is zero when k > K. Other filters can 
be used, such as 

= («) 

where e is recognised as a regularisation parameter, and the jk are eigenvalues of a 
constraint operator. These are the basics of solutions via SVD. In this solution no 
assumption has been made about the operator H except that it is positive definite. 
In the case of a spatially invariant blur this operator has what is known as a Block 
Toepliiz structure. This can be further approximated so that the structure becomes 
block circulani i.e. the' operator becomes periodic, and solutions to the imaging 
problem can then be found via FFTs. The solution of most imaging problems are 
found by making these approximations. The SVD method does not, and therefore 
allows the inversion of operators formed with a spatially variant PSF. As we have 
already mentioned conventional SVD upon reasonably sized images is at present 
impossible. The scanning method allows a route to the solution of this problem. 


5.3 The scanning method 

To overcome the speed and virtual memory requirements of SVD upon large images 
we have made use of the fact that in many problems the PSF is very local in nature. 
We also make the assumption that the PSF has a finite support within the image 
space. Then a far smaller portion of the singular system is required to reconstruct 
the image. If an appropriately sized portion of the image is taken it can be used to 
reconstruct the corresponding part of the object via conventional SVD. 

g, = H.f, + T], (49) 

In this equation the image vector g, corresponds to a small portion of the full image 
and Hi, fi and T]i are the corresponding operator, object portion an noise vector. 


32 



In equation 49 the operator will be far smaller e.g. 1024 X 1024 for a 32 X 32 
pixel object and image region, which allows an SVD in a reasonable time. In fact 
the SVD on the operator for the 32 x 32 image region will be O (10 ) times quicker 
than an SVD on an operator for an image size of 512 X 512. After the singular 
system is obtained the reconstruction via equation 49 is obtained in the order of 
n'* operations. The reconstruction for the 512 pixel image will require of the order 
512“' operations whereas a 512 pixel image reconstructed in 32 x 32 pixel parts will 
require of the order of 32“ X 16^ operations. This gives a factor of O (10^) in speed 
up. Therefore there is approximately a speed up of 0 (10’°). 

To reconstruct the image in parts as suggested above is rather crude and the 
reconstructed object will contain errors due to what we call blocking effects due 
to the overlap of the PSF in neighbouring reconstruction areas. To overcome this 
problem we simply reconstruct in a smaller box so that image points outside the 
image box have no effect. The scanning box is the portion of the image used, the 
reconstruction box in chosen to lie central in the scanning box. To reconstruct the 
whole image the scanning box is shifted so that the next reconstruction box lies 
next to the previous one. In this manner the blocking effects are avoided and good 
reconstructions can be obtained. All that is lost is an area around the edge of 
the image which cannot be reconstructed as information is needed from outside the 
image. It must also be noted that the factor of 10^° speed increase wiU be reduced 
by this procedure. 

As an improvement on this we also introduced an averaging method which in¬ 
stead of reconstructing in small boxes next to one another the reconstruction boxes 
were allowed to overlap. An average of the reconstructions was then taken. This 
allowed improvements in reconstructions as numerical errors were partially removed 
by the averaging process. 

The operator for the segmented part of the image in equation 49 still has had 
no further assumptions made about its form. Therefore the spatially variant PSF 
problem can stiU be tackled with the scanning technique. In the case of a spatially 
invariant system only a single SVD has to be performed and the singular system 
can then be stored and,'used again. In the spatially variant PSF many SVDs have 
to be performed which means that the method becomes computational rather slow, 
but it is still accurate. In fact within this smaller area of the image the PSF should 
have a weakly spatially variant nature. A possible solution is then to assume spatial 
invariance within this area and use Fourier transforms to solve which wiU be far 
more computationaUy efficient. 

It may also be noted that equation 49 also shows that if the image has a spa- 
tiaUy variant noise upon it, and that noise can be considered as constant across 
the scanning box, then treatment of this problem is quite straightforward. AU that 
is required some knowledge of the signal to noise ratio across the image so that a 
regularisation procedure can be suitably varied when reconstructing each portion 
of the image. This was not investigated further within this project but would be 
interesting to foUow up 

5.4 Results and discussion 

To show the blocking effects mentioned in the previous section a scanning SVD was 
performed upon the simple image of a cross convolved with a Lorentzian of radius 
one pixel. Figure 30 shows the results and the errors in the reconstruction are given 


33 



in table 2. In figure 30a where the reconstruction is performed with scanning and 
reconstruction boxes of the same size the blocking effects are obvious and give a high 
reconstruction error. As the size of the reconstruction box is reduced the blocking 
effects disappear, and a good reconstruction is obtained with a low error as shown 


in figure 30f. 



(a) Reconstruction box of size 16 x 16 pixels, 43.0% error, (b) Reconstruction box 

of size 14 X 14 pixels, 6.0% error. 



(c) Reconstruction box of size 12 x 12 pixels, 5.0% error, (d) Reconstruction box 

of size 10 X 10 pixels, 5.0% error. 



(e) Reconstruction box of size 6x6 pixels, 3.0% error, (f) Reconstruction box of 

size 2x2 pixels, 0.5% error. 

Figure 30. The variation of reconstruction error with reconstruction box size. 
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Table 2: Errors in reconstruction with reconstruction box size 

To show the performance upon a more complex image a 512 X 512 image of Lena 
figure 31a was obtained and convolved with a Lorentzian of radius one pixel. The 
image has 5.0% Gaussian white noise was added. This image can be seen in figure 
31b. A reconstruction was obtained using the scanning method with a scanning box 
size of 32 x 32 pixels and a reconstruction box size of 2 x 2 pixels, figure 31c. In this 
case 157 of the available 1024 singular values were used in the reconstruction. 





Figure 31b. Image with Lorentzian PSF of radius 1.0 pixels, condition number 

138, and 5.0% Gaussian noise. 



Figure 31c. A reconstruction 30.0% in error with the true image. The 
reconstruction was performed with a scanning box of size 32 x 32 pixels, and a 
reconstruction box of 2 x 2 pixels with 157 singular values. 
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The real advantage of this method is that it can be used to reconstruct images 
that are formed with a generally spatially variant blur. As mentioned earlier to do 
this with SVD will be very time consuming, and an approximate method which uses 
the same scanning idea can be performed with Fourier transforms. In this case it is 
then assumed that across the scanning box area that the PSF is spatially invariant. 
This is of course an approximation. An experiment was performed with this method 
upon an image blurred by a Lorentzian whose width increases with the distance from 
the centre of the image. The results are shown in figure 32. 



Figure 32a. The original object. 



Figure 32b. An image formed with a spatially varying PSF. The PSF is Lorentzian 
in form with radius 1.0 pixels at the center increasing linearly to 3.0 pixels at the 

edge. 
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Figure 32c. The reconstructed image using the approximate Fourier scanning 

method. 

The image shown is noiseless and when initial reconstructions were performed 
large errors occurred. This 'was due to the locally spatially in'variant PSF approx¬ 
imation enforced by the use of Fourier transforms. To cope with this problem a 
the regularisation parameter was given a non-zero value. The departure from the 
spatially invariant case is then treated as noise, and controUed in the same manner. 
The use of SVD in reconstructing the image with spatially variant blur gives good 
results although it is rather computationally slow. In this case the regularisation 
parameter does not need altering to compensate for departures from the local spatial 
invariance of the PSF. An example reconstruction is shown in figure 33. 


5.5 Conclusions 

We have introduced a scanning method for SVD reconstruction of images. This 
method removes the blocking effects that would normally be associated with a re¬ 
construction by parts. It is also vastly superior in terms of computation speed as 
compared with conventional SVD. The method is also applicable to the case of spa¬ 
tially varying blurs. The results have shown that the blocking effects are removed 
to a tolerable level with the scanning technique, and good reconstructions can be 
obtained with large complex images containing noise. To show the potential for 
the reconstruction of images formed with a spatially variant PSF an approximate 
method using Fourier transforms was shown. This introduced artifacts into the re¬ 
construction which were controlled by use of a regularisation parameter. An example 
reconstruction has been shown for an SVD reconstruction which is more accurate 
than the Fourier method, but far slower in computational speed. 




Figure 33a. The image of Lena with a spatially variant blur. 



Figure 33b. The reconstruction. 









6 Further work 

In this contract we have investigated a number of methods for obtaining diffraction 
limited imaging when the optical systems are degraded, for application to such 
problems as for example space to space imaging. We have applied three lines of 
research to this problem, (i) Blind deconvolution, (ii) Image reconstruction using a 
known aberration function, and (iii) Exponential-filter speckle imaging. Our major 
effort has been upon the first two of these and we have made significant contributions 
to these fields. 

We have investigated the performance of existing blind-deconvolution algorithms; 
particular issues have been robustness of algorithms to noise and the object/PSF 
morphology. We have also investigated a number of novel approaches to blind- 
deconvolution. (a) Blind deconvolution using the Richardson-Lucy (Expectation 
maximisation of maximum likelihood (EMML)) algorithm, and (b) semi-blind de- 
convolution, which makes use of further a priori information, i.e. having a known 
form for the PSF but unknown parameters describing it [10]. 

In the blind-deconvolution work some significant results have been obtained with 
the use of both full and semi-blind Richardson-Lucy algorithms [10]. In the blind 
case, good reconstructions have been obtained with as much as 15.0% photon noise in 
the image. In the semi-blind case a PSF was chosen with three free variables, where 
the functional form of the PSF was chosen to be a Gaussian plus a Gaussian times 
the radius squared. The free parameters were then the central Gaussian width, the 
ring width and the ring height. These parameters must be found by the algorithm. 
A non-linear least squares fitting routine was used at each step to find the PSF 
parameters, and it was found that the algorithm had convergent properties with as 
much as 4.0% noise upon the image. This is a great improvement upon previous 
blind-deconvolution work. 

Image reconstruction with a known aberration Tu.nction has been mainly per¬ 
formed upon HST images before the corrective optics were added. In terms of the 
BMD project the HST provided real data which would have similar deficiencies to 
any other orbiting optical system. A number of algorithms were considered i.e. 
Richardson-Lucy, conjugate gradients and Landweber with and without and pos¬ 
itivity constraints, and Tikhonov regularisation [13]. The large size of the HST 
images made it impossible to use conventional singular value decomposition (SVD) 
methods. To overcome this limitation we have recently introduced a new technique 
which scans the image and produces reconstructions using SVD. The scanning SVD 
method reduces the computation time for image restoration by orders of magnitude 
without detriment to the quality of the image by exploiting the linear and^ local 
nature of the PSF. The technique is potentially of great importance because it can 
be used for both isoplanatic and non-isoplanatic imaging This will allow further 
improvement of the HST images even with corrective optics, and has application to 
a whole host of other imaging problems. 

Experimental work has been done in the speckle imaging area. A photon lim¬ 
ited imaging system with variable turbulence degradation has been built. It has 
been found in the past that detectors have been unable to cope with moving speckle 
patterns when trying to form the auto-correlation of the pattern due to ‘dead time 
effects. A detector with good ‘dead time’ characteristics was loaned to us by Ruther¬ 
ford Appleton labs (RAL). It was found that at suitable light levels the detector 
gave good auto-correlations with an acceptable level of noise for the moving speckle 
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pattern. It is now possible therefore to obtain real data in the laboratory to test 
algorithms. 

To summarise in the current contract we have made a significant breakthrough 
with the new scanning SVD techniques. We have made important contributions 
to the technique of blind-deconvolution and we have made useful comparisons of 
performance for existing blind and non-blind deconvolution algorithms. We have 
also set up a satisfactory test bench for generating data. 

Areas in which we envisage that this work could be continued is based upon 
the progress made with the SVD techniques for efficient image restoration in high- 
numerical-aperture systems, and to continue the work on blind deconvolution with 
the use of a priori information. Again these subjects have direct applicability to high 
resolution space imaging as well as for microscopy, astronomy and medical imaging. 

For general high-numerical-aperture imaging our current work on scanning SVD 
techniques is directly applicable. In the case of isoplanatic imaging the scanning SVD 
technique has made possible for the first time the use of SVD on large images, and 
has made a large increase in the speed at which modest images can be reconstructed. 
Due to the scanning nature of the algorithm it can be generalized to high-numerical- 
aperture systems where the isoplanatic condition no longer exists. 

One advantage of using SVD is that the singular system once calculated can 
be stored and high-resolution imaging performed for any object at speed. The 
scanning technique further improves upon this. In the non-isoplanatic case further 
improvements are still required. It is proposed that theoretical and computational 
investigations be made into the efficient SVD of matrices which are similar but not 
the same, which is the case when performing a scanning SVD upon an image formed 
by a non-isoplanatic system. 

Other methods for the efficient reconstruction of images formed by optical sys¬ 
tems with a spatially variant blur would also be of interest. One particular area 
which would be relevant to the scanning SVD techniques we have developed, is the 
concept of ‘Displacement Rank’ [14]. This allows the efficient inversion of Toeplitz 
matrices, which are thq operators formed by a one dimensional spatially invariant 
blur. The Toeplitz structure can be relaxed and these methods still work, there¬ 
fore a weakly space variant PSF system can also be inverted efficiently. The use of 
the scanning technique and this method should allow more general spatially variant 
systems to be tackled. 

Another area of interest is the preconditioning of imaging operators for fast 
convergence. In the work of Plemmons [15], for example circulant preconditioners 
have been used upon block-Toeplitz operators [16], i.e the operator formed by a 
two-dimensional spatially invariant system. The preconditioning effectively clusters 
the singular values of the operator around unity so that algorithms such as conju¬ 
gate gradients converge very rapidly. The investigation of preconditioning for the 
spatially variant PSF imaging system complement our work on SVD. 

In the blind deconvolution area further investigations into the amount of a priori 
knowledge required for convergent algorithms would be a good area for further inves¬ 
tigation. In all the algorithms to date except the semi-blind algorithm demonstrated 
in the current work no guaranteed convergent algorithms have been produced. The 
latter algorithm incorporated a large amount of a priori information in the func¬ 
tional form of the PSF. Research could be carried out to investigate the amount 
of a priori information required to give a convergent and reliable algorithm. The 
deconvolution process is reliable in the sense that the algorithms are convergent and 
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give results that, depending upon the amount of noise present, give a good repre¬ 
sentation of the object. If the problem is altered so that a priori information is 
available for both the object and PSF, how much information is required to give a 
reliable algorithm. To perform this some factorisation of probability distributions 
is required. Linnik has made significant contributions to this field [17,18] and it is 
hoped that these methods could be applied to the above problem. 

Finally it may be practical to combine the fast linear scanning SVD method with 
the non-linear blind-deconvolution or other methods by use of the output values of 
the former as starting data for the latter. Advantages would be the satisfaction of 
constraints which can built into the non-lineaj methods without the large number 
of iterations normally necessary. 
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